subroutine w_recompute_FM_mags_step_1(iter,c,nd,stanum,evnum,nsta,nev, &
           p,maxdist_P3,fm_corrs,inverted_mags,   &
           MLFM_sta,MLFM_final_res,ML_ev,ML_sd,sigma2)
! 


    implicit none
    include 'inversione_pesata_multipar.FI'
    include 'seisev_EVENT_PHASE_structure.FI'
   
    integer km,ks,k1,k2,ir,ic,ier,ndimcur,iter
    integer k,nd,nsta,nev,i,j,kk,idgt,ie
    real*4  inverted_mags(ievents),fm_corrs(istations),rms
    real*4  bstar(npartot),wkarea(nnw),resdato,somma
    integer stanum(150000),evnum(150000)
    real*4  MLFM_sta(150000),ML_ev(5000),ML_sd(5000)
    real*4  MLFM_s_wgt(150000),MLFM_final_res(150000)
    real*4  maxdist_P3,p(6)
    real*4  PAI,MLFM1,MLFM2
    real*8  sigma2
     real*4 prova,ssdd(5000),nndd(5000) ! deviaz standard e numero di fasi
   
    integer last,id,first
    logical flag
    integer*8 memo
    real*4  v(5000),wmag,wsd,wmax,wmin,ww(5000),ruse
    real*4  ML, A, R, CS, MLFM, res, whuber, MAXDIF_EVSTA
!   A  = ampiezza WA              memorizzata in c(i).WA_pha_amp
!   R  = distanza ipocentrale     memorizzata in c(i).s_h_hypdist
!   CS = correzione di stazione   memorizzata in c(i).DBcorr

!   Dichiarazione di function:

!   MLDB(A,R,CS) = log10(A) + 1.667 * log10(R/100) + 0.001736 * (R-100) + 3 + CS

!    MLFM(A,R,CS) = log10(A) + P1 * log10(R/100) + P2 * (R-100) + 3 + CS &
!    + P3 *SIN ( (((R  /100)**2.)**0.31)  *2.0 *PAI )  * &
!           EXP(-((R  /200)**0.5) ) 

!    + P3 * sin(SQRT(R)*pai/5.) * exp(-R/600.)

    MLFM1(A,R,CS) = log10(A) + P(1) * log10(R/100) + P(2) * (R-100) + 3 + CS &
    + P(3) * sin(  ((R/100)**P(4))  *2.*pai)         * exp(-(R/P(5))**P(6))

    MLFM2(A,R,CS) = log10(A) + P(1) * log10(R/100.) + P(2) * (R-100.) + 3. + CS

    PAI=3.1415927
    
    write(*,*)' subroutine: recompute_FM_mags'
    !   Write magnitudes out:
    open(unit=6,status='unknown',file='FM_event_magnitudes.csv') 
    write(6,'(a)')'ot,lat,lon,depth,mag,FM_ML,sd,n-dati,ML_Rich,ML_DB'   

    
!   percorro tutti i dati 
!   e calcolo per tutte le stazioni una magnitudo FM di staz.
!   Ho delle correzioni per tutte le stazioni.
    do i=1,nd
      if(c(i).WA_pha_amp .le. maxdist_P3)then
        MLFM_sta(i)= MLFM1(c(i).WA_pha_amp,c(i).s_h_hypdist,fm_corrs(stanum(i)))
      else
        MLFM_sta(i)= MLFM2(c(i).WA_pha_amp,c(i).s_h_hypdist,fm_corrs(stanum(i)))
      endif
    enddo


!    do i=1,nd
!        MLFM_sta(i)= MLFM(c(i).WA_pha_amp,c(i).s_h_hypdist,fm_corrs(stanum(i)))
!    enddo
    
!   Ora tocca calcolare la magnitudo d'evento.
    first = 1
    last  = 1
    do while( evnum(last) .eq. evnum(first) )
      last = last +1
    enddo
    last = last -1
    kk = 0
    ie=0
    do while (last .le. nd  .and. first .le. nd)
      ie = ie +1
      k= 0
      do i=first,last
          k=k+1 
          v(k) = MLFM_sta(i)
      enddo

      ruse = 0.23
      wmag = 0.
      wsd = 0.
      call Huber_dev(k,v,ruse,wmag,wsd,wmax,wmin,ww)
            kk=kk+1
      ssdd(kk)=wsd
      nndd(kk)=float(k)

      ! memorizza i risultati   !
      !

      
      flag = .FALSE.
      k=0
      ML_ev(ie)=wmag
      ML_sd(ie)=wsd
      do i=first,last
          if(evnum(i) .ne. ie)then
             write(*,*)' errore in event number. i, evnum(i), ie:',i,evnum(i),ie
             pause
          endif
          k=k+1
  !        MLFM_ev(i) = wmag                            !event mag
  !        MLFM_sd(i)   = wsd                           !event maf standard-dev
          MLFM_final_res(i) = MLFM_sta(i) - ML_ev(evnum(i))  ! residuo dopo inversone
          MLFM_s_wgt(i) = ww(k) ! whuber(res,flag,ruse)   ! peso della mag di stazione    
          c(i).MLRH_s_wgt = ww(k) ! whuber(res,flag,ruse) ! ridefinisco peso per nuova invers.   
      enddo   
           
      ! OCCORRE RICALCOLARE I PESI ASSEGNATI A CIASCUNA ML DI STAZIONE PER
      ! EVENTUALMENTE RIUTILIZZARLI DI SEGUITO NEL CALCOLARE LE CORREZIONI DI STA.
      ! RIPERCORRO I DATI. IN QUESTO PASSO RICALCOLO I PESI PER TUTTE LE STAZIONI,
      ! SIA PER LE STAZIONI DIBONA CHE PER LE ALTRE
      !  flag = .FALSE.
      !  ruse = 0.3
      !  do i=first,last
      !     res = c(i).MLDB_sta - c(i).MLDB_ev
      !     c(i).MLDB_s_wgt  = whuber(res,flag,ruse)
      !  enddo
       
        j=first
        write(6,'(1x,a16,a1,3(f7.3,a1),3(f6.2,a1),i5,2(a1,f6.2) )') &
          c(j).ot,',',c(j).lat,',',c(j).lon,',',c(j).depth,',',  &
          c(j).mag,',',wmag,',',wsd,',',k,',',c(j).MLRH_ev,',',c(j).MLDB_ev
        
     
      
      ! vai al prossimo evento
      first = last + 1
      last = first
      do while( last .le. nd .and.  evnum(last) .eq. evnum(first))
        last = last +1
      enddo
      last = last -1

    enddo        
    close(6)
 ! rewrite all the data with new outcomes
    open(file='FM_mags.csv',status='unknown',unit=7)
 
    write(7,'(a)')'sta,cha,ot,s_h_hypdist,amp1,amp2,wa_pha_amp,'//   &
                   'MLDB_dif,MLDB_sta,MLDB_s_wgt,MLDB_ev,MLDB_sd,'// &
                   'lat,lon,depth, mag, N_WA_DB,NEV,N_WA_P_EV,'//    &
                   'N_WA_P_STA,stanum,evnum,inverted_mags,ML_ev,'// & 
                   'ML_sd,MLFM_sta,MLFM_s_wgt,fm_corrs,MLFM_final_res,MLRH_sta'    
 
    sigma2 = 0.0d00
    do i=1,nd
      sigma2 = sigma2 + dble(MLFM_final_res(i) * MLFM_final_res(i)) 
      write(7,'(a5,a1,a3,a1,a19,a1,4(f15.5,a1),9(f7.3,a1),6(i6,a1),       &
                8(f7.3,a1))')                                             &
      c(i).sta,',',c(i).cha,',',c(i).ot,',',c(i).s_h_hypdist,',',         &
      c(i).amp1,',',c(i).amp2,',',c(i).wa_pha_amp,',',c(i).MLDB_dif,',',  &
      c(i).MLDB_sta,',',c(i).MLDB_s_wgt,',',c(i).MLDB_ev,',',c(i).MLDB_sd,',', &
      c(i).lat,',',c(i).lon,',',c(i).depth,',',c(i).mag,',',c(i).N_WA_DB, &
      ',',c(i).NEV,',',c(i).N_WA_P_EV,',',c(i).N_WA_P_STA,                &
      ',',stanum(i),',',evnum(i),',',                                     &
      inverted_mags(evnum(i)),',',ML_ev(evnum(i)),',',ML_sd(evnum(i)),',',  &
      MLFM_sta(i), ',',MLFM_s_wgt(i),',',    &
      fm_corrs(stanum(i)),',',MLFM_final_res(i),',',c(i).MLRH_sta   
    enddo
    write(*,*)' sigma2 ',sigma2
 
      rms=SNGL(sigma2)
      rms=rms/float(nd)
      rms=SQRT(rms)
      
      
      write(11, '(a,a,i3,a,f10.2,a,i5,a,f8.4,a,f5.3,a,f9.6,a,f9.6)') &
 'My inversion 3 par',',',iter,',',sigma2,',',nd,',',rms,',',P(1),',',P(2),',3.0,',P(3)

      
      
    do i=1,kk
      prova = prova + ssdd(i)*ssdd(i)*nndd(i)
    enddo
    
    WRITE(*,*)' = = = = = P R O V A = = = = = MIA STEP 1  ',prova
    

    
 return
 end     
 
 
 
 